install.packages("dplyr")
install.packages("DESeq2")
install.packages("ggplot2")
install.packages("clusterProfiler")
install.packages("org.Dr.eg.db")
install.packages("enrichplot")
install.packages("DESeq2")
install.packages("stringr")
library(dplyr)  
library(DESeq2)  
library(ggplot2)  
library(clusterProfiler)  
library(org.Dr.eg.db)  
library(enrichplot)
library(DESeq2)
library(stringr)
data=read.csv("C:/Users/Rommel.Dks/Desktop/r/GSE217196_GEOcounts.txt", header=TRUE, row.names=1)
histdata <- data
hist(log2(as.matrix(histdata)),main = paste("对数化直方图"))
library(dendextend)
library(ggtree)

#3
data1 <- t(data)
data1 <- scale(data1)
data2 <- dist(data1)
data3 <- hclust(data2)
plot(data3)
#4
EOF_list=grep('^(EOf)',colnames(data))#8个样本
MPS_list=grep('^(MPS)',colnames(data))#8个样本
wt_list=grep('^(wt)',colnames(data))#8个样本
EOf_DEG <- cbind(data[EOF_list],data[wt_list])
EOf_DEG <- EOf_DEG[which(rowSums(EOf_DEG)>10),]#过滤掉EOf组中低表达值基因
MPS_DEG <- cbind(data[MPS_list],data[wt_list])
MPS_DEG <- MPS_DEG[which(rowSums(MPS_DEG)>10),]#过滤掉MPS组中低表达值基因
EOf_group <- ifelse(substr(colnames(EOf_DEG),1,3)=="EOf","EOf","wt")
EOf_group <- factor(EOf_group,c("EOf","wt"))
MPS_group <- ifelse(substr(colnames(MPS_DEG),1,3)=="MPS","MPS","wt")
MPS_group <- factor(MPS_group,c("MPS","wt"))
DEG <- EOf_DEG
group <- EOf_group
coldata <- data.frame(row.names =colnames(DEG),
                      condition=group)
dds<-DESeqDataSetFromMatrix(countData = DEG,colData = coldata,design = ~condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition","EOf","wt"))
res <- data.frame(res)
write.csv(res, file = "分析结果.csv")
res <- res[complete.cases(res),]#去除含有缺失值的行
res1 <- subset(res, pvalue < 0.05 & abs(log2FoldChange) > 1.0)
write.csv(res1, file = "筛选分析结果.csv")
#火山图
data4<-read.csv("C:/Users/Rommel.Dks/Desktop/r/分析结果.csv")
colnames(data1)[1] <- "ID"
data4$log10pvalue <- -log10(data4$pvalue)
data4$significant <- ifelse(data4$pvalue < 0.05 & data4$log2FoldChange > log2(1.2), "Up",
                            ifelse(data4$pvalue < 0.05 & data4$log2FoldChange < log2(1/1.2), "Down", "Stable")) ###设置上调、下调和稳定的条件
#绘制火山图
p<-ggplot(data4, aes(x = log2FoldChange, y = log10pvalue, color = significant)) +
  geom_point(size=1) +
  scale_color_manual(values = c("Up" = "red", "Down" = "green", "Stable" = "black")) +
  theme_minimal() +
  xlab("Log2 Fold Change") +
  ylab("-Log10 P-value") +
  ggtitle("EOf_Volcano Plot") +
  theme(legend.position = "right")
ggsave("EOf_volcano_plot.jpg")
plot(p)
#go
data5<-read.csv("C:/Users/Rommel.Dks/Desktop/r/筛选分析结果.csv")
colnames(data5)[1] <- "ID"
data5$ID[1:5]
gene_id<-data5$ID
entrez_genes <- bitr(gene_id, fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Dr.eg.db)
go_enrich <- enrichGO(gene = entrez_genes$ENTREZID,
                      OrgDb = org.Dr.eg.db,
                      keyType = "ENTREZID",
                      ont = "ALL",
                      pAdjustMethod = "BH",
                      pvalueCutoff = 0.05,
                      qvalueCutoff = 0.05)
head(go_enrich)
dotplot(go_enrich, showCategory = 20) +
  ggtitle("EOf_GO Enrichment Analysis") +
  theme_minimal()
ggsave("EOfdotplot1.png", width = 10, height = 8)
barplot(go_enrich, showCategory = 20) +
  ggtitle("EOf_GO Enrichment Analysis") +
  theme_minimal()
ggsave("EOfbarplot1.png", width = 10, height = 8)
#7分析
DEG <- MPS_DEG
group <- MPS_group
coldata <- data.frame(row.names =colnames(DEG),
                      condition=group)
dds<-DESeqDataSetFromMatrix(countData = DEG,colData = coldata,design = ~condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition","MPS","wt"))
res <- data.frame(res)
write.csv(res, file = "MPS分析结果.csv")
res <- res[complete.cases(res),]
res1 <- subset(res, pvalue < 0.05 & abs(log2FoldChange) > 1.0)
write.csv(res1, file = "MPS筛选分析结果.csv")
data7<-read.csv("MPS分析结果.csv")
colnames(data7)[1] <- "ID"
data7$log10pvalue <- -log10(data7$pvalue)
data7$significant <- ifelse(data7$pvalue < 0.05 & data7$log2FoldChange > log2(1.2), "Up",
                            ifelse(data7$pvalue < 0.05 & data7$log2FoldChange < log2(1/1.2), "Down", "Stable")) 
p<-ggplot(data7, aes(x = log2FoldChange, y = log10pvalue, color = significant)) +
  geom_point(size=1) +
  scale_color_manual(values = c("Up" = "red", "Down" = "green", "Stable" = "black")) +
  theme_minimal() +
  xlab("Log2 Fold Change") +
  ylab("-Log10 P-value") +
  ggtitle("MPS_Volcano Plot") +
  theme(legend.position = "right")
ggsave("MPS_volcano_plot.jpg", plot = p, width = 8, height = 6, dpi = 300)
plot(p)


